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20  ABSTRACT  (Continued) 

frequencies  above  i'  10  Hz,  and  by  practical  computer-time  constraints  to  frequencies 
belr\v  v -100  Hz,  The  code  allows  the  calculation  of  the  amplitude  and  phase  of  the 
full  acoustic  field  as  a function  of  depth  and  ran^o  at  any  desired  time. 
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I INTRODUCTION 


Sonar  has  I, eon  used  for  detection  and  communication  in  the  ocean 

for  many  years.  The  effect  or.  sonar  of  a geographical  variation  in  the 

so, ,„d -speed  profile  of  the  ocean  has  received  considerable  attention  as 

more  precise  and  sensitive  use  has  been  made  of  acoustic  signals  in  the 

ocean,  however,  relatively  little  attention  has  been  focused  on  the 

smaller- scale  variations  in  the  ocoans-in  particular,  the  internal 

waves,  which  not  only  provide  variations  in  space  (on  scales  of  a few 

hundred  meters  vertically  and  a few  kilometers  horizontally),  but  also 

introduce  the  new  dimension  of  time  variability  <o„  a scale  of  a few 
hours ) . 

The  theory  of  sound  propagation  through  internal  waves  has  not  yet 
been  fully  developed,  so  that  many  question,  of  a qualitative  nature  are 
open,  we  have  written  a computer  code  to  answer  some  of  these  questions, 
taking  advantage  of  the  advent  of  the  use  of  the  parabolic  equation  method 
acoustics,1"1*  and  the  existence  of  internal-wave  models 
derived  from  oceanographic  measurements. “ With  certain  restrictions  (to 
be  described)  the  full  acoustic-wave  field  can  be  propagated  through  an 

ocean  with  variable  sound  speed,  and  the  resulting  long-range  transmission 
characteristics  can  be  studied. 


References  are  listed  at  the  end  of  the  report. 


Con. si dci  the  acoustic  pressure  in  the  ocean  as  p(r,t).  The  wave 
equation  Cor  p is 


2 

V P 


2 

-1 

2 2 
c at 


0 


where  it  has  been  assumed  that  the  soiind  speed  c(r,t)  varies  slowly  with 
position  and  that  the  water  density,  as  far  as  its  effect  on  inertia  is 
concerned,  is  essentially  constant.  Hence,  the  quantity  that  determines 
the  structure  and  behavior  of  the  pressure  is  the  sound  speed  c(r,t). 

In  the  following  sections  we  present  a model  for  the  expected  sound- 
speed  structure  in  the  ocean,  the  parabolic  equation  approximation  to  the 
wave  equation  and  the  reasons  it  is  valid  for  our  model  of  the  sound-speed 
structure,  and  the  computer  code  that  has  combined  the  model  of  internal 
waves  and  the  parabolic  equation  method  into  a useful  apparatus  for  studying 
oceanic  sound  transmission. 

Initial  results  obtained  through  use  of  our  code  are  given  in  Ref- 
erence 6. 
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II  OCEAN  SOUND-SPEED  STRUCTURE 


A ■ Dotorminis tip  Profile 

On  the  scale  of  the  depth  of  the  ocean  (4  to  5 km)  the  sound  speed 
as  a function  of  depth,  z,  is  determined  by  the  gross  behavior  of  the 
density,  temperature,  and  salinity.  We  use  the  profile  derived  by  Munk,5 
whose  input  is  an  exponentially  decreasing  density  gradient.  The  result- 
ing "canonical"  profile  is 


c (z) 
CPV  ' 


= 1 + 6 [e_11  - (1-7))]  | 


where  Tl  = —r~ 


Note  that  c(z)  has  a minimum  at  z , that  the  width  of  the  minimum  is  B, 

and  that  the  deviation  of  the  sound  speed  from  the  minimum  value  c is 

1 

of  order  e . 

Typical  (though  not  universal)  values  for  the  parameters  are: 

I 

zA  = 1300  m 
A 

B = 1300  m 


e = 0. 7d  x 10 


e = 1500  m/s 


See  Figure  1 for  a graph  of  ccp(z)  using  these  parameters. 


Illustrations  arc  grouped  at  the  end  of  the  report. 
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One  must  point  out  that  realistic  ocean  profiles  in  most  cases  have 
significantly  different  behavior  from  this  general  form.  For  example, 
within  a few  hundred  meters  of  the  surface  a mixed  layer  usually  results 
in  a lowered  sound-speed  gradient.  However,  we  ignore  these  details  in 
the  present  version  of  the  code. 


B,  Internal  Waves 

The  density  gradient  in  the  ocean  leads  to  the  possibility  of  waves 
traversing  the  volume  of  the  ocean  just  as  the  density  discontinuity  at 
the  surface  leads  to  tine  possibility  of  surface  waves.  The  density 
gradient  is  usually  presented  in  the  form 


N (z) 


( Po6z  ) 


where  N(z)  is  called  the  local  stability  (Brunt-Vaisala)  frequency. 

i 

Munk's  model  for  N(z)  is 


N ( z) 


-z/B 

n e 
o 


where  n =4  cycles/hr. 
o 

Let  w(r,t)  be  the  displacement  of  a particular  density  level  from 
its  mean  position.  It  can  be  shown  that  w satisfies  the  equation: 


The  eigenmodes  of  this  equation  can  be  found  by  taking 

iTkjX+kgy-w ( j ,k)t] 

w = W ( j , k , z ) e 
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k W = 0 


Modifying  this  equation  to  account  for  the  rotation  of  the  earth  we  fi.i 


u>  -to 


k W = 0 


where  u>  = inertial  frequency  = (2  cycles/day )sin (latitude) . 

Boundary  conditions  arc  W(z)  = 0 at  surface  and  bottom.  Thus,  for 

a given  horizontal  wave  number  there  are  a discrete  infinity  of  solutions 

to  this  equation  labeled  by  the  mode  number  j,  with  frequencies  on ( j , k) . 

The  sound-speed  fluctuation  profiles  caused  by  each  internal  wave  will 

2 

follow  the  depth  structure  of  N (z)W(j ,k,z) . Several  examples  of  such 
profiles  are  shown  in  Figure  2. 

The  sound  velocity  fluctuations  caused  by  a full  internal-wave  field 
will  be 


Nw  = Ro 


Uk1x+k2y-<i)(j,k)t]  „ 
k2)o  Nz(z)W(j,k,z) 


J’kl'k2 


Ihe  difficulty  of  integi’ating  this  distribution  over  the  x-z  plane 
(because  we  will  propagate  sound  only  in  the  x-z  plane  vith  the  acoustic 
transmission  code)  has  caused  us  to  consider  a simplified  version  of  the 
internal-wave  spectrum  where  internal  waves  are  propagated  only  in  (or 
opposite  to)  the  direction  that  the  sound  waves  propagate.  In  addition, 


. - " -V*  f ^ ■ t '■  S 


wc  combine  real  and  imaginary  parts  to  reduce  fluctuations  in  the  overall 
energy  in  the  internal  waves  as  a function  of  time: 


'1W  yfe 


He  A r 1 in  A 


A = ^ A(j,k)  W(j,k,z)N2(z)  eA^kr 

j,k 


where  W(i 


'1,k,z)  is  »ormalizod  so  that  J m“V  (z)W2(j  ,k.z)  dz  = 1. 


The  A ( j , k ) are  complex  Gaussian  random  variables.  From  a synopsis 
of  diverse  oceanographic  measurements,  Garrett  and  Munk4  have  proposed 
the  following  model: 


<A«. kWh' 


4cu  k^j  (cy-1)  (j  . f'1 


o / flUU.  j 

, . . ,cv  , 2 I i 

n (J+J*)  k + - 

o * \ n B 

' o 


4 x 10 


whoie  = J,  Qf  = 2.5,  and  L is  the  horizontal  periodicity  length  over  which 
the  internal  waves  are  generated  (because  we  are  generating  k in  discrete 

steps  rather  than  as  a continuum).  The  relative  intensity  of  the  various 
modes  is  shown  in  Figure  3, 


* tv. : «.  * 7, 
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N'olo  t.liat.  the  magnitude  of  the  sound-speed  fluctuation  due  to  in- 
ternal waves  is  6c/c  ~ 10  , a factor  ol  one  hundred  be  low  t lie  magnitude 

of  the  deterministic  structure.  Also  the  spatial  behavior  of  the  sound- 
speed  variations  due  to  internal  waves  is  of  the  order  of  a few  hundred 
meters  vet  tically  and  several  kilometers  hori/.ontall'  . Internal-wave 
frequencies  vary  from  n = I eyeles/hr  for  the  lowest  mode  to  one  eyelo/day 
for  the  high  modes  (,j  ~ 20).  Figure  1 shows  some  typical  sound-speed 
profiles  due  to  internal  waves. 


C . Final  Expression  for  Sound-Speed  Structure 

The  sound  speed  as  a function  of  position  and  time  is  given  by 
o(T-, t)  = c <z)  *tclt 


where 


0(10“4) 


Hence,  the  deviation  in  sound  speed  from  a constant  value  is  always  small, 
and  we  may  write,  for  example, 


This  expression  will  of  course  appear  in  the  wave  equation. 
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III  THE  1MKAB01.IC-EQUAT10N  METHOD 

A . Int  rod uc  L i ? .11 

l h<  pa  ''a  bo]  ie-equa  t.ion  method  was  originally  developed  by  Ixjontovich 

and  I- ok  in  19dG  Lo  study  long-range  propagation  of  radio  waves  in  the 

troposphere.  This  method  was  introduced  into  the  field  of  underwater 

acoustics  by  Tappcrt  in  1972,  and  a computer  program  based  on  this  method 

uas  developed  by  Tappcrt  and  Hardin  to  solve  acoustic  propagation  problems 

2 3 

of  interest  to  the  Navy.  ’ 

•I.  Approximations  and  Ranges  of  Validity 

1 ho  wave  equation  for  acoustic  pressure  p(r,t)  is 

r,-  1 3~p 

p ■ T 3 = 0 • 

c d t 

Since  the  time  variation  of  c takes  place  over  hours,  and  the  acoustic 
frequencies  will  be  in  the  range  of  ~ 200  Hz,  we  may  treat  sound  propa- 
gation at  any  time  as  though  the  ocean  structure  were  frozen. 

donee,  a particular  acoustic  frequency  uu  would  be  unaffected  by  the 
ocean  and  we  may  sot 

, -ioi  t 
P = P e 

2,2,  a) 

V p + k p =0  where  k 


9 


but  k 


2 


2 , 2 , 2 6c  , 

- V p i-  k p -2k  — p = 0 . 

o o c 

Notice  tint;  the  c/c  term  initially,  the  solution  in  cylindrical  co- 
ordinates is 


! k,, )' 


p oc  ||  (k  r)  — — 
o o \r 


so  we  try  for  the  full  solution: 

p'  = \J/  (r,/,.*p)ll  (k  r) 
t o o 

where  the  reduced  wave  function  \J/  is  labeled  by  the  time  t,  because  the 
c c structure  of  the  ocean  is  different  for  different  times.  Substitute 
int;  p'  in  the  full  equation  we  find: 

o ’>  5ll 

II  + W”ll  + 2 ~ — + k2II  V - 2k2  — II  Y = 0 

o o or  or  o o o c o 


If  k r » l, 
o 


dll  / 

° / 1 

= i k 1 1 + 0 

or  o l 2ik  r 2 

\ o k 

O 


b)  "» 


• V2Y  - i — + 2ik  2 *£  t . o 

r o r o o r o c 


.2  2 2 

3 V 1 3 t a V 


3V 


c 


, + —IT  + — - J 2ik  — - 2k  — Y = 0 
• ~ 2 2 2 o dr  o c 

o r r sep  a z. 


The  key  to  the  parabolic-equation  method  involves  the  following 
physical  approximations,  based  on  the  structure  of  tjc/c: 


10 


• k r » 1 
o 


2 i k 

o 


y± 

Sr 


\alid  if  the  objects  off  which  the  acoustic  waves  are  scattering  have 

sizes  that  are  much  larger  than  a wavelength.  This  is  equivalent  to 

having  only  relatively  forward  scattering,  which  results  in  small  changes 

in  \|'  over  an  acoustic  wavelength.  Lot  L be  the  vertical  scale  of  sound 

v 

speed  variations.  Then  k^L^  » 1 is  required.  Note  that  energy  scattered 
at  largo  angles  is  removed  by  the  ocean  bottom  anyway. 


2 2 2 
• LI  _L  LL  - L_i 
2 2 2 ' 2 
or  r Scd  Sz 

Tills  is  true  because  the  canonical  profile,  which  is  100  times  larger 
than  the  sound-speed  fluctuation  of  the  internal  waves,  affects  the  z 
coordinate  only.  More  importantly,  however,  the  internal-wave  gradients 
in  the  vertical  arc  an  order  of  magnitude  greater  than  the  horizontal. 

The  approximate  wave  equation  is  therefore 


1 


2k 

o ciz 


i 0 c 

+ k — \|/  = 

o c 


(1) 


To  summarize  the  approximations  required  for  this  parabolic  equation  to 
be  valid,  we  need  the  following  quantities: 


(i) IW  - Largest  frequency  involved  in  the  internal-wave  spectrum 
53  4 cycle  s/hr 

L = Minimum  horizontal  scale  of  sound-speed  fluctuations 


Ly  - Minimum  vertical  scale  of  sound-speed  fluctuations 
«=  200  m. 

Validity  of  the  parabolic  eciuation  requires: 

( 1 ) 0)  » U) 

I IV 

(2)  k r > 1 

o 

(3)  L » L 

11  V 

(d)  k L » 1 . 
o V 

ltequi roment  3 is  well  satisfied  by  our  internal-wa ve  model.  Re- 
quirement 1 sets  a lower  limit  on  the  frequency  we  can  treat: 

1 C1 

ireq  >->  — — « 1 II  z . 

V 

Requirement  1 is  subservient  to  Requirement  -1.  Requirement  2 sets  a 
lower  limit  on  the  range  of  sound  transmission: 

R » — 25  in  for  10  Hz 

k 

o 

Thus,  any  experiments  with  frequencies  of  10  Hz  or  higher  over  ranges 
of  (10)/k(j  or  longer  should  be  well  treated  by  the  parabolic  equation. 

This  includes  all  problems  of  interest  at  present. 


Analogy  with  Quantum  Mechanics 

The  Schrodinger  equation  in  nonrclativistic  quantum  mechanics  is 


i 

- + V(z)t 

2m  2 


ass 


Bt  ' 


(2) 


12 


We  see  that  our  parabolic  equation  tl)  is  of  the  same  form  ii  the 


analogies  are  made  as  follows; 

Mass  = k 

o 

f-c 

Potential  = k — 
o c 

Time  = range 

Energy  level  = horizui'al  wave  number 
/ dz 

Momentum  = angle  — * 

\ llr 

Many  characteristics  of  underwater  acoustics  can  be  easily  understood 
on  the  basis  of  this  analogy.  We  will  give  two  examples  for  transmission 
with  6c /c  given  by  the  canonical  profile  (no  internal  waves): 

(1)  Rays — Consider  a narrow  wave  packet  with  velocity  v propagating 

2 

in  our  potential  well  T see  Figure  5(b)!  . As  long  as  |?mv  < V 
the  packet  will  bounce  back  and  forth  inside  the  well.  The 
analogy  says  that  as  long  as  -g-k09  « kn  (6c/c)max,  a ray  start- 

r 

ing  from  t ho  sound  axis  with  0 will  be  contained  within  the 
sound  channel  (see  Figure  5(a)’]. 

(2)  Ground  Ktate--Forming  a wave  packet  with  the  form  'j/Q(z)  where 

\j'0  is  the  ground-state  wave  function  of  the  potential  V(z)  would 
result  in  an  acoustic  field  propagating  in  range  with  no  change 
in  structure — i.e.,  no  spreading  of  the  wave  packet  (see  Figure  6). 
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1\  NUMERICAL  REALIZATION  OF  THE  INTERNAL-WAVE  MODEL 


A . Eigenfunction  Determination 

7 

A program  (/MODE)  was  used  to  numerically  generate  the  eigenfunctions 
W(j,k,z)  and  frequencies  o>(.j,k).  Several  modifications  to  ZMODE  were  re- 
tin  i red : 


(1)  An  increase  of  dimensions  so  that  128  values  of  k and  256  values 
of  7.  can  be  treated. 

(2)  Inclusion  of  the  inertial  frequency  oi  . 

i 

(3>  Oblige  of  the  initial  guess  Joru>(.1,k).  |y^ (0)  = instead 

of  (m-A ) it . I 

c I 


B.  Range  of  ,j  and  k 


Modes  with  1 < .j  •?  2d  are  included.  Values  of  k range  from  (k  /128) 

max 

to  (127  I28)kmax  in  127  equal  steps.  Both  positive  and  negative  values  of 

k are  included,  making  a total  of  25'1  values  of  k.  'Hie  value  of  k is 

max 

1 cycle/km  at  present.  Thus  a total  of  254  X 24,  or  6096  separate  modes, 
are  generated. 


C.  Generation  of  6c 

1W 


Ihc  steps  in  generating  the  soun.i-speed  deviations  caused  by  internal 
waves,  6c^,  are  as  follows: 


(1)  Ihc  6096  different  A(j,k)  arc  generated  according  to  a 
Rayleigh  probability  distribution  in  amplitude,  and 
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variance  Riven  by  the  sped  rum  described  in  Section  1J.  The 
phase  angle  oj  cacli  A(j,k)  is  randomly  generated  in  the  region 
0 t°  - • Thus,  12192  random  numbers  arc  needed. 

lor  each  ,j  value,  the  values  ol  k are  stepped  through,  calcula- 
tion of  the  function  W(.j,k,z)  is  accomplished,  and  a matrix 
clement  is  incremented  until  the  following  complex  matrix  is 
completed: 


M ( k , z ) = EA(j,k)S“(z)IV(j,k,z)c1*(j,k)t 

j 

The  value  of  6c  at  range  r and  depth  z is  the  Fourier  transform 

of  M(k,z):  u(r,z)  = S M( k , z)e , The  matrix  A is  calculated 

1 rcim  M by  the  last -Fourier^t  ransform  technique.  The  values  of 

r at  which  A is  known  arc  0 < r <-  256  /k  in  stens  of  -/V 

max  * max1 

In  the  present  code  this  means  0 < r 128  km;  the  code  simply 
repeats  the  structure  it  has,  every  128  km. 

We  repeat  that 

c j 

5°,,,,  = — [RoA  + ImA]  . 

1 U 2 

Hence  we  have  the  internal-wave  sound-speed  deviations  at 
time  t on  a grid  of  128  km  by  4 km  with  256  points  horizontally 
and  256  points  vertically.  (Step  size  ~ 16  m vertically  and 
0.5  km  horizontally.) 

The  propagation  section  of  the  code  may  need  the  value  of  6c 

IW 

on  points  off  the  above  grid.  The  code  allows  division  of  the 
grid  into  a finer  grid  by  powers  of  two,  with  the  values  of 

6oIW  t,n  the  rino  grid  being  obtained  by  cubic  spline  interpola- 
tion (see  Appendix), 


(5)  The  25(5  256  grid  of  bc^,  values  is  kept  in  Large  Core  memory 

ior  fast  access  during  the  sound-propaga  t ion  phase  of  the  code. 
When  the  values  of  6cjw  at  a new  time  t are  desired,  the  entire 
procedure  is  repeated  with  a new  lime  t,  with  two  exceptions; 
first,  the  random  numbers  do  not  need  to  be  generated  again, 
and  second,  the  eigenvalues  u,(j,k)  do  not  need  to  be  calculated 
again  because  they  have  been  saved  in  computer  memory.  Due  to 
storage  limitations,  however,  the  functions  W(,j,;;,z)  must  be 
recalculated  at  every  t,  since  they  have  not  been  saved. 


V NUMERICAL  PROPAGATION  OF  THE  ACOUSTIC  FIELD 
USING  THE  PARABOLIC  EQUATION 


A.  Algorithm 


Wo  must  numerically  solve  the  parabolic  equation 
^ f 1 d"  6c  w 

57  = 7T  - lk„  T * = 1<A  + B,Y 

o dz 


wheio  reference  to  the  variable  z is  suppressed,  and  where  A and  B are 
real  operators.  The  first-order  solution,  useful  for  a numerical  method, 


dr  dr 

, i x iA  -n-  iBdriA-o- 
V (i+di)  = e ^ e e 2 \J/(i-) 


Since  step  follows  step,  wc  use 

m,  , , , iAdr  iBdr.„ 

* (r  I-  dr)  = e e Y (r) 

.....  iBdr 

-valuation  of  e is  straightforward — it  is  a complex  number  with  unit 
modulus.  Evaluation  of  c is  more  complicated  since  A is  a differen- 
tial operator: 


iAdr  n ^ iAdr 

c ri(z)  = F e Flti(z)  ] 


F is  a fa s t - Four ic i'- 1 ransf orm  opei’ation,  and 


— — k2 
2ko 


in  k-space, 
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This  algorithm  is  fast,  and  also  it  is  very  stable,  sine 


since  the  total 


aeons t i c one rgy 


/W2"z 


>s  exactly  conserved  as  a lunction  of  r.  This  is  because  the  operations 
<>n  ( r)  are  exactly  unitary. 


h.  Acoustic  Source 

The  acoustic  field  may  be  started  with  any  function  of  depth.  The 

present  code  models  a point  source,  at  a depth  provided  as  input,  with 
unit  strength  at  1 yard. 

C . Surface  Boundary  Condition 

The  reduced  wave  function  vj/  must  be  exactly  zero  at  the  ocean  sur- 
face (pressure  release  boundary  condition).  In  order  to  ensure  the 
proper  application  of  this  condition,  a "reflected  ocean"  is  carried 
along  with  the  real  ocean.  In  the  reflected  ocean  is  an  image  source 
of  the  same  strength  as  the  real  source,  but  reflected  in  Position  and 
180°  out  of  phase  with  the  real  source  (see  Figure  7). 

f* • Bottom  Absorption 

In  order  to  avoid  reflection  off  the  ocean  bottom,  a gradual  loss 
of  intensity  is  imposed  on  *(z>  as  z nears  the  ocean  bottom.  The  func- 
tional form  of  the  imposed  loss  is 

< ...  ,I,.\  , , v iAdr  iBdr 

V ( i +di  ) = h(z)e  e \J/  ( r) 


lAz)  = exp  { - a dr 


The  code  has  n t.  present 


Ot  0.05/m 

e = o.o4  /. 

max 

This  form  effectively  stops  any  acoustic  energy  from  penetrating 
below  about  500  in  above  the  bottom.  Even  this  form  does  cause  some 
reflection  off  the  bottom  at  a low-intensity  level,  for  long-wavelength 
sound . 

E . Number  of  Depth  Points 

The  maximum  number  of  depth  points  in  the  real  ocean  allowed  by  the 
code  is  1024  (the  number  used  must  be  a power  of  two). 

F.  Size  of  Range  Step 

The  size  of  the  range  step  should  not  be  too  large  or  the  numerical 
algorithm  will  break  down.  Experience  has  indicated  that  0.5  km  is  a safe 
step  at  100  Hz  and  0.25  km  is  safe  at  400  Hz. 


\ 
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VI  FORMAT  OF  SAVED  DA'I’A  FROM  COMPUTER  RUNS 

As  t lie  sound-propnga  t ion  code  is  run,  various  quantities  are  put 
on  magnetic  tape,  much  in  t ho  manner  of  an  actual  field  experiment.  The 
tapes  can  then  he  separately  analyzed. 

At  the  present  'he  output  consists  of  the  following  logical  records: 

(1)  An  identification  record  (see  Table  1). 

(2)  A number  of  data  records  containing  the  acoustic  field  at 
various  points  (see  Table  2).  For  each  time  step  the  full 
acoustic  lield  is  output  lor  12  values  of  the  range , and 
the  ncousric  field  at  four  selected  depths  is  output  for 
all  values  of  the  range. 

(3)  An  end-of-filc  mark. 


r 


Table  1 


IDENTIFICATION  RECORD 


Word  Number 


Quan 1 i t v 


Forma t 


Table  1 (Concluded) 


Word  Number 


Qunn lily 

Forma t 

NF,  (internal  waves) 

I 

niinx 

. N K (internal  wn  ve  s ) 

1 

I 

NMODE  (internal  waves) 

1 

1 <(6c/c)2)  | z=0] 

F 

T I IV 

F 

•i* 

I 

Number  of  range  steps 

1 

Number  of  time  steps 

I 

Time-step  interval 

F 

Table  2 
DATA  RECORD 
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VII  COMPUTER  RUNNING  TIMES 

The  breakdown  <>!'  the  CDC  7600  compute!’  time  to  run  the  eodc  is  as 
f ol lows : 

(1)  To  generate  the  6c  on  a 256  X 256  grid  for  each  time: 

I\V  

28  CPU  seconds. 

(2)  To  step  the  acoustic  field  once  in  range 

0.13(N  log  N)'/(10  X 1024)  CPU  seconds  where  N is  the 
number  of  deptli  points  at  which  the  acoustic  field  is 
ova lua ted. 

Hence  the  total  running  time  is 

-5 

T = NT 1 28  + 1.3  X 10  (N  logi;)N)NRJ  CPU  seconds 

where 

NT  = Number  of  time  steps 
N = Number  of  depth  points 
NR  = Number  of  range  steps  at  each  time. 

Note:  As  of  September  1974  the  code  has  been  changed  to  make  use  of  a 

fast  sine  transform,  and  the  eigenfunctions  W(j,k,z)  are  being  saved  on 
disk.  This  has  resulted  in  a faetor-of-two  increase  in  speed. 
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YIN  CODE  ORGAN!'/  AT  ION  AND  LISTING 


Figure  8 illustrates  the  program  flow  and  t lie  FORTRAN  listing  of 
the  code  in  its  present  form  is  given  in  the  following  pages. 
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PROGRAM  USERIOUTPUT  , T A PE6 =GUTPUT , T APE l 0 > 

COMMON  C ( 30  I 
0 1 MENS  I ON  I C ( 301 
EQUIVALENCE (C,  I c I 
DIMENSION  OPR ( 8 ) ,DBT( 8 » 

COMMON/NEACTS/N,LN,NR,NTLGSS,NIT 

CGMMON/AFACTS/FREG,SNORM,ZMAX,GW»FGW,DR»FKR,DZ»WL,Pl»GD 

CGMMON/VAR/R, IR,T, I T 

COMMON /AC OUST / UR ( 20481, Ul(20 48) 

COMMON/ STUP/ISTOP 

COMMON/ PH A SES/SK (1024) ,CK| 102^) ,SM( 102^1 ,CM( 102^1 
COMMON/ BFUR/UR 81 1 024 * , U I 8 ( 1 024 > , ANGB(  1024) 

COMMON/ MUNK/C A, 2 A, ZW, GEE 

COMMON/ ImAVES/ENO, FINER, OZ,EZMAX,NZ,AKMAX,NK, NNQOE  »DELC  * 

* TIw.JSTAR 

COMMUN/ SAV/SAVRI 100*4) ,SAVl I 100,41 
COMMON/ STEP /NS TEP, I ST E P ,NSAV ,DT 
DIMENSION  LR(12),LZ(4) 

DATA  LR/20, 40, 60, 80, 100, 120, 160, 200, 2t>0, 320, 400, 5 00/ 

DATA  LZ/25, 102, 166,256/ 

DATA  ENO,FlNER,NZ,AKMAX,NK,NMODE,DELC,T iw,jstar/ 

* 4, ,u. 04, 257, 1.0, 128, 24, 0.0, 2. 5, 6/ 

DATA  T, IT/O., 1/ 

DATA  CA,ZA,ZW,GEE/1500.,1300.,1300.,0.017/ 

DATA  ISTOP/O/ 

DATA  FRED  *GD, DPR , S NORM/ 100. ,1300.  ,1000. ,7*0. ,1500./ 

DATA  ^gpk^^NR,NTLUSS,ZHAX’FGw  /5  12,  .25 ,200, 1 ,4000.  , 1 ./ 

DATA  R , RKM  /2*0./ 

DATA  PI  /3. 1415926356/ 

R = 0 . 

Z BF  = 1 JOO . 

I R=  1 

GU  TO  60 
70  CONTINUE 

CALL  T APE  I W C 1 1 
DO  50  I T = 1 , 1ST  EP 
MR=  1 
NSA V=0 
R = 0. 

NP  = 0 

T =D T* ( IT-1  I 
CALL  INIT 
CALL  PHASE1 
CALL  PHASE  2 
DO  10  1R=1,NSTEP 
JK=MOO( IR.mOI 
IF  (UELC.EQ.O. I GO  TO  767 
CALL  PHASE2 
767  CONTINUE 
CALL  UWS 

r=k«-or 

IF  (MODI IR.201 .EU.01  CALL  PR  I NT P( UR , U I , N , K » 

NSAV=NS A V ♦ 1 
DO  14  IDE=l,4 
I DP  =LZ ( I DE  I 

SAVR(NSAV,IDE)=UR(  I DP  I 
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U SAVI(NSAV,1DEI=UI(  I 0 P I 

IF  UR  .NE  .LRIMR  » I GO  TO  12 
MR=MRf 1 

CALL  TAPE  I w ( 2 ) 

NS AV=0 
12  CONTINUE 
TO  CONTINUE 
50  CONTINUE 

WR  1 T fc ( fa , 2 00  I 
I S T CP  = 1 

CALL  TAPE IW ( 3 ) 

STOP 

GO  CONTINUE 
NS  T E P=  500 

FRFU=100.  i N = 5 12 
DPKH-O. 5 
OT  = 3840. 

N S A V = 0 
OE  L C= . 0004 
I S TEP=  32 

BZ  = ZW  > E ZMAX  = ZMAX 
T P I = 2 . *P  I 
EN0=ENQ*TPI/3600. 
FlNER=FlNEP*TPI/3faQ0. 
AKMAX=AKMAX*TPI /1000. 

W8F  =2  50 • 


LN-ALOGIFLOAT ( N I I/AL0GI2. ) ♦ . 5 
N=2**LN 

OR=ORKM* 1 000.  i WL=SNORM/FREQ  $ FKR=2.*Pl/WL 
Gw=FGW*WL  i UZ = ZM A X /F L OAT ( N I 
ZKM=ZMAX/ 1000. 

RANGE  = NS  tep*drkm 

*WKDELC 6,3001  ,N,DZ,RANGE,NSTEP,DRKM 

300  FORMAT! lHl,*  ACOUSTIC  FREQUENCY  *,F4.0,*  HZ*// 

* OCEAN  DEPTH  *,F3.l,*  KM  NO.  OF  DEPTH  POINTS*. 15. 

* STEP  SIZE*, F 8. 2.*  M*// 

* RANGE  *,F6.2,*  KM  NO.  OF  RANGE  STEPS*, 15, 

* STEP  SIZE*,F5.2 ,*  KM*// 

* UELC  * , E I 2 . 2, / / ) 

WR I TE ( 6, 100 i 

WR I TE ( 6 , 200 ) 

100  FORMAT ( 20X, *PLOT  OF  THE  ACOUSTIC  FIELD*// 

* 17X,*ALL  DIMENSIONS  ARE  IN  KILOMETERS*,// 

* * RANGE  / DcPTH  ->*/ 

* * ?*/ 

* * ?*) 


200  FORMAT! 1H0,9X»*0.0*9X*0,5*9X*I.Q*9X*1.5*9X*2.0* 
* 9X*2.5*9X*3-0*9X*3.5*9X*4.0*/I 
GO  TO  70 
END 
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SUBROUTINE  E lot  NFNl  GAMMA,  K,W,WG,0Z  ,NZ,FZN> 
COMMON  /EGNFN/  F ( 2 57 I , G( 2 57 I , ON ( 2 57 » , UNS ( 2 57 1 
REAL  K f KU 

DATA  C6,C2A/.l  6606666b  66  b 66  , • 0 A 166  6660600  66  / 
PU(X> =GAMMAG*X-KU 
GAMMA J= (DZ*GAMMA)**2 
KU= (DZ*K  I **2 
F ( n=Fl=FZZI  = 0. 

FZ1=0Z 

G(  l MG1=GZ1=GZZI  = 0. 

00  luo  1=2, NZ 
F2  = FU.5*FZU.125*FZZ1 
G2  = G 1 ♦ • 5*GZ 1 ♦ • 12  5*GZZ 1 
P02  = PQ  ( UN  ( I-lll 
F2=F2-C2A*(FZZ1*F2*PU2> 

U2  = G2-C2A*1GZZI«-G2*PU2»F2*UN{  I- II  ) 

FZZ2=-F2*PU2 

GZZ2=-G2*P02-F2*0N l I-l I 

F3=F1*FZI*C6*(FZZ1»2.*FZZ2I 

G3  = G1«-GZI*C6*(GZZ1 *2.*GZZ2> 

PU'J  = PO(  ONI  I > I 
FZZ  3=-F3*PJ3 
GZZ3=-G3*PU3-F3*UN( I I 
F ( I I =F 1 = F 3 
G<  I I =G 1 = G3 

FZl=FZl*Cb*(FZZl*4»*FZZ2*FZZ3l 

GZI=GZ1+C6*(GZZUA.*GZZ2*GZZ3) 

FZZ1=FZZ3 
100  GZZ1=GZZ3 
W*-F  1 

*G=-Gl*0Z**2 

FZN=FZ1 

RETURN 

END 
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SUbROU I I NL  £ I GENVL I GAMMA , K , NORM , JNORM , MODE , MERR , I ERR, 
.NZ,DZ,HN, LRR) 

COMMON  /FONFN/  F { 2 57 ) , G ( 2 5 7 » , ON < 2 57 I , QNS ( 2 5 7 ) 

COMMON/ 1WAVES/ENO.F INER,AI9) 

REAL  K,  NORM 
MERR=0 
I E RR  = U 
DGMAX= .o/HN 
SGN=l . 

I F ( MOO{MOOE,2 I .EQ.OI  SGN=-l. 

DO  10  1=1,10 

CALL  E IGENFNIGAMMA,K,W,WG,DZ,NZ,FZN) 

I F ( SGN**G.LE.0.1  GQ  IQ  11 
DGAMMA=W/(2.*GAMMA*fcG ) 

IF  ( A*3  S < UGAMMA)  .GT  . 0OMAX  ) DGAMMA  = S ! GNI  DGMAX,  DGAMMA  1 
GAMMA=GAMMA-OGAMMA 

IF { AdS I UGAMMA/GAMMA I .LE  .ERR  I GO  TO  12 

10  CUNTINUF 

11  I ERR= 1 

12  M=  1 

EPS=0Z**2*W/,<G 
DO  20  I =1 ,NZ 
20  Fin--F||)-EPS*G(I) 

N1  = NZ-  1 
OU  25  1=2, N1 

IF{  (F(  I-1)*F'U  » .LI  .0.1  .UR.  Fill  .EQ.O.I  M = M*1 

25  CONTINUE 

MERR  =M-MOl)E 

IFIMERR.Nfc.O.OR.IERR.EO.l » RETURN 

ONURM=~rtG*FZN/DZ 

RON=l./SUkT (ONURM) 

SUM=0 . 

SUMR  = 0 • 

DU  30  1 = 1, NZ 

SUM=SUM*F( I ) **2 

F(  I )=F  ( I | *(  UNI  II  ♦ F INER**2  I 

SUMR=  SUMR  *■  F(  I 1**2 

30  CONTINUE 

UNURM= 1 . / SURT ( SUMR*DZ) 

DO  31  l=l,NZ 

31  FI  I ) = F(  I )*UNURM 
NORM=UZ*SUM 
RETURN 

END 


SUttkUU  TINE  FFr842lN2PUk,X,Y) 

DI  MENS  I UN  X(2l,Yt2l,HI5l,CS(2G4  9>,SS(2  0A9> 

EUUIVALENCE(Ll5,L(in,(LlH,L(2»»,(L13,L(3H,(L12*L(41J,<Lll,L(5M, 

ULl0,Ltfall,tL9,Lt7IML8,L(8M,(L7tL(qil,a6,Lll0lh<L5,Ulll).IL<t 

2«L(  12  1 1 , ( L 3 , L ( lit  ) * ( L 2 *L  ( 1 A 1 I , ( L 1 , L ( 15)  ) 

DATA  LSI  /O/ 

NTFIPU  =2**N2POW 
N8Pr)H  = N2POW/  3 
I F ( N8P0W  ) 17*3.17 
17  LNG=NTHPU 
FLNG=LNG 

IF  (LNG-LSU  171,173,171 
I 7 I L S T = LNG 

DT=6. 283185 30 7179/FL MG 
CSC  11  = 1. 

SS( 11=0. 

TH  = 0. 

MD=LST/8 
DU  172  1=2, MD 

ih=th*dt 

csi  n=  cost  TH» 

172  SSI1I=  SIN(TH) 

173  DO  1 I PASS® 1 , N8P0W 
NXTLT  =2**(N2POW-3*IPASSI 
LENGT  =8*NXTL  T 
SC=FLMG/FLUAT(LENGT  ) 

D01 J=1 .NXTLT 

MR  G = F L UA I ( J-l I * SO ♦ 1 .5 

C 1 — C S C MRG  1 

S 1 =S  S ( NRG  I 

C2=C1**2-S1**2 

S2  = Cl*Sl«-Ll*Sl 

C3=Cl*C2-Sl*S2 

S3=C2*S  1 + S2*C 1 

CA=C2**2-S2**2 

SA=C2*S2*C2*S2 

C5=C2*C3-S2*S3 

S5  = C 3*S2«-S3*C2 

Co=C3**2-S3**2 

So  = C3*S3«-C3*S3 

C7=C3*Ca-:>3*SA 

S7  = C4*S3«-SA*C3 

DOIISULO  =L  ENGT  , NTHPQ  ,L  ENG  T 

JO® I SOLO  -LENGT  *J 

Jl=JO>NXrir 

J2  = Jl*’ NXTLT 

J3* J2* NXTLT 

JA= J3+NXTLT 

J5= JA+NXTLT 

Jti=  J 5*  NX  TL  T 

J7=J6+ NXTLT 

AKU=X( JO  I +X( JA  I 

AR 1 = X ( J 1 ) ♦ X( J 5 * 

AR2=X< J2I+XI Jo» 

AR3*X( J 3 » ♦X C J7) 

ARA  = X( J 0 » - X ( JA » 

AR5  = X< J1 » - X < J5) 
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ARb  =X ( J2  » - X ( JbJ 
AR7=X( J 3 » - X < J 7 I 
A I 0 = Y ( J01  *-Y(  J4> 

AII=YU1)*Y(J5) 

AI2  = Y(J2»4-Y(J6I 
AI3  = YCJ3»«-Y(J7* 

AU--YI  JO»-Y(  J*| 

A I b = Y ( J 1 I - Y ( J b » 

Al6=Y( J 2 I - Y C J6) 

A l 7 = YC  J 3 I - Y I J 7 I 
BRO=ARO* AR2 
BR  1 = AR  UAR3 
flR2=AR0-AR2 
BK 3 = AR 1 - AR  3 
BR^»  = AR4- A I 6 
BR8  = ARt>-A  l 7 
BRo=ARA>AI6 
BK7=AR5*Al 7 
B10  = AI0*-AI2 
B I 1 = AI  I *A  13 
B I 2- A 1 0- A I 2 
B I 3 = A I 1-AI3 
B I 4 = A I <»♦■  AR6 
dI5=AI5*AR7 
B16=AI 4-AR6 
B I 7 = A I 5-AR7 
X(  JOI  =BR0*BR1 

Y ( JO  I = B I 0*B  It 
IFU-ll  2,2,18 

x ( J 1 1 = C4*(  BR0-BR1>-S^*JRI0-Bin 
Y(  J1)=C4*(  BIO-BI  i)*$4*IBR0-BRU 
X( J2 I =C2  * C BR2-BI3»-S2*(BI2*BR3) 
Y( J 2 I = C2 * ( BI2»BRJ|»S2*(BR2-8I 31 
X(  J3I  = C6*(  BK2*BI3|-S6MB12-BR3> 
Y(J3)=Cb*(BI2-BR3l*S6*(dR2*Bl3l 
TR=  0.7071067812*1 BR5-BI5) 
TI=0.7071067812*(8R5*BI5I 
X(J<M=Cl*(BR<H-TRl-Si*CBU«-m 

Y (J4l  = Ci*(  B l4*Tn*Sl*(BR<r*rR| 

X(  J5I=C8*  ( BR^-TR|-S5*<  BH-T  I ) 
Y(J5>=C8*(BI4-TI t ♦! 5*( BR^-TR  1 
TR=-0.7071007812*(BR7+BI7 ) 
11=0.707106781 2* (bP 7- B17> 
X(J6)=C3*(BR6*TRI-S3*<BIbHH 
Y(J6|=C3*(BI6*T! I ♦ S3  * C BR6+TRJ 
X(J7»=C7*(BR6-TRI-S7*( BI6-TI) 

Y( J7 1 =C  7* ( B lo- 1 1 I ♦ S7* ( BRb-TR I 
GO  TO  1 

X< J1)=BR0-BK1 

Y ( J 1 1 =8  I 0-  B I 1 
X( J21=BR2-B  I 3 
Y ( J2>=B12*BK3 
X(  J3  1 = BR 2*B  l«3 
Y( J 3 1 = B I 2-BR  3 

TR  =0.707  1 0670  12*1 13R5- B I £»  I 
T 1=0.  70  71  06  7812*(  8R5*-B  181 
X(  J4)  = BR<**TR 


Y(  J9  > = Bl<»»Ti 
X( J6»=BR9-TR 
Y( J6 I =8 l 9-T  I 

TR--0.70710to7612*(BR7*Bl7  > 
TI=0.7071067B12*(UR7-BI7» 
X( Jb )=BR6+TR 
Y( Jbl=B16*T  I 
X( J7I=BR6-TR 
Y( J7»=BI6-Tt 
1 CUNTtNUE 

3 IF(N2POW-3*N8POW-1I 6,6,7 

6 00b  1 J=  1 , NTHPO  ,2 

R1=X(JI*X(J*1I 
X( J»1I=X( J I -X  4 J ♦ 1 1 
X( JI=R1 

Fll=Y(JI*Y(J*l) 

Y( J»1)=Y(J)-Y( J ♦ 1 > 

61  Y I J I =F  1 1 

GOT  U5 

7 UU7 1 Jl  = l *NTHPO  ,4 

J2>JU1 
J 3 = J 2 ♦ 1 
J4-*  J 3*  1 

R1=X<  J 1 I 4-X  C J3  I 
R2*X(J1»-X( J3> 

RJxX(  J2I  *X(  J'*) 

R4=X(J2 l-X( J9I 
FI1*Y(J1)*Y(J3I 
F 1 2 = V C J1J-YU3) 

F l 3 = Y ( J2 I ♦ Y ( J4  J 
FI4=Y(  J2»-Y( J4> 

X(  JH=RUR3 
YIJ1)=F11*F13 
X( J2I-R1-R3 
Y( J2J  = F Il-F  13 
X( J3I=R2-FI9 
Y ( J 3 I = F 1 2 +R4 
X(  J<»I=R2*F  l<» 

71  Y ( J<*  I =F  I 2-R^ 

5 DU6 1 J= 1 » 15 

L ( J ) = 1 

I F ( J-N2PUW ) 19,19,61 
19  L ( J I = 2**( N2POW* 1-J I 
61  CUNTINUE 
I J*  1 

DUB  J 1= 1 , L 1 
D08J2=J 1 ,L2,L 1 
DUB J3= J2 , L3,L2 
DU8 J9= J 3, L9 ,L  3 
D08J5  = J<*,L5,L^ 

008  J6  = J5,L6,L5 

DD8 J 7 = J6 , L 7 , L6 

0U8J8=J7,L8,L7 

DQ8J9=J8,L9,L8 

D08J10=J9,L10,L9 

D08J11=J10,L11 ,110 

D08J12«=J11,L12,L11 


OOdJ13=Jl2,Ll3,L12 
008  J14  = J 13  ,L  14  ,L  1 3 
0U8JI=J14,L18,L14 

I F ( l J- J | | 20,8,8 
20  R=X(IJ) 

A(  I J)  = X( Jl  | 

X ( J I ) = R 
F I = Y ( I J ) 

Y(  I J»=Y(  Jl  ) 

Y(  JI)=M 
8 I J * I J 4- 1 

RE  TURN 
END 


SUBROUTINE  I N I T 

COMMON  / AC'OUST  / UR  (2048)  , U I (2048) 

COMMON  /NFACrS/  N,LN,NR,NTLQSS,NIT 

S3??!:,  “ 

S G A I N x ( .9144/Gw) *SORT( ML/PI  ) 
ou  10  J=1,N 
QEPTh=DZ*FLOAT (J) 


UR(  J)=PRUE(L)EPTH)-PROF(-DEPTH) 
10  UI ( J)=0. 

RETURN 

ENO 


10 


SU8RUU TINE  PHASE  1 

cSJmS  /ofattq/  ^ 1024), CM  1 024  ) , SM(  1024), CM  ( 1024) 

CUMMUn  /nfacts/  n,ln,nr,ntloss,ni t 

SNOHM*  ZMAX»GW»FGW»0R«FKR»D2,yL»Pl»G5 

FN=FLUAT ( N*2 ) 

00  10  K= 1 , N 


AK=FLOAT ( K ) *FK/FKR 

ANG=.5*DR*FKR*AK**2 

SK(K)=SIn(ANG)/FN 

CK(K)-COS(ANG)/FN 

RETURN 

ENO 
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SUBRUU T I Nfc  PHASE2 

REAL  SSM(  102<i  ) * SSF  ( 1024ltALOSSl  102A> 

s;: 

DATA  ISTRT  /O/ 

I F ( ISTRT  .GT.  01  GOTO  100 

I S T RT= 1 


DLW=.OA*ZMAX 
DL  P= . 0 5 


DO  10  J=1,N 
DEPTH=DZ*FLOAT ( J < 

BLOSS3DLP*EXP( -( ( DtPTH-ZMAX ) / DL  W ) **2l 

10  ALOSSC  J)=EXP(-0R*BL0SS) 

DO  20  J=1,N 
OEPTH=DZ* FLOAT ( J ) 

ARG  = 2.*<  OEPTH-ZAI/ZM 

20  SSM(J)*CA>GEE*(ZR/2. )*(ARG~1.+EXP(-ARG)I 
DO  JO  J=1,N 
EN=SN0RM/SSM1 J I 

30  SSM( J1=.5*FKR*DR*( EN»i.|*(EN-l . > 

100  CONTINUE 

CALL  RNDMF ( SSF 1 
DO  120  J = 1 » N 
ANG=SSM( J I ♦SSF { J I 
SN( J1=AL0SSIJ)*SIN( ANG) 

120  CM( J)=ALOSSC J) *CUS( ANG) 

RETURN 

END 


SUBROUTINE  PRINTP 

COMMON/ AC 0 US T/UR( 20  48 ) ,UI(20A8> 

COMMON/ NF  ACTS/NtLN,NRtNTLOSS»NIT 

CCJMMUN/VAR/R,  IR 

INTEGER  P ( 32  I 

1 SK=N/ 32 

RAN=R/1000. 

DO  10  J=l,32 
I“( J-l 1*1 SK*l 
P0W=UR< | 1 **2*U I ( I ) **2 
10  PC  Jl=-l0.*AL0G10(P0W/Rfl.E-100> 
W«ITEC6,20l  RAN,P,IR 
20  FURMAT1F7. 2, 3X, 3213,151 
RETURN  $ END 
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FUNCTION  PS(K,MQDEI 
REAL  K 

COMMON  / I W A VE  S / ENU»FINER|6Z»LZMAX»NZ»AKMAX,NK,NM0DE» 
I S IG, T IW, JSTAR 
OATA  ISTRT  /O/ 

DATA  PI  IS. 1415S26336/ 

H(X»=(T1*-1.I*(1.+X»**(-TIW> 

IF ( ISTRT. Nk. 01  GO  10  100 
I S TR  T = 1 
S=  F INfcR/ENO 
AKO*l./BZ 

AA=FLOAI (JSTAR I *P  I *S*AKO 
AK1=AKMAX/FL0AT(NKJ 
CQEF  = AKI*S1 0**  2 * ( BZ  / 3 • I *(4.*S)*AK0 
100  ALAM=FL OAT (MODE  I /FLOAT (JSTAR! 

AKRFF=AA*AL AM 

PS=CUEF*K**2*ALAM*M(  AL AMI /( K**2*AKREF**2 1**2 

RETURN 

END 
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SUBROUTINE  RNDMF(SSF) 

REAL  SSFt It , SSTI257) ,YZ(257) ,AZ1257) 

COMMON  / 1 WAVES/  ENO . F I NER , BZ . EZMAX , NZ . AKMAX , NK , NMODE . 

I S IGt  T I M i JSTAR 

LARGE  EFR(256,257) , EF  I ( 2S6,  2 57  ) 

COMMON  /NFACTS/  N , L N , NR ,N T LOS S , N I T 

COMMON  /AFACTS/  FREQ , SNQRM, ZMAX , GW ,FGW , DR ,FKR , DZ , WL , P I 
COMMON  /VAR/  R, IR,T,  IT 
REAL  ESK(256),ESI (256) 

REAL  SP{4) 

DATA  ITT  /O/ 

I.F(  IT.EU.ITH  GOTO  IOO 
I TT  * I T 
DO  1 I =1  ,N 

1 SSF (I  ) = 0 . 

IFISIG.EU.O. ) GOTO  IOO 
NX=2*NK 

TP  1 = 2 . *P I 
□ X=  P I /AKMAX 

SR2»1./SQRT(2. » s 

FCTR*FKR*DR 

UEZ  = EZMAX/FLOAT(NZ-I  ) 

R0fcZS=l./DEZ**2 
CALL  Z MODE J 

LNX=ALUG (FLOAT ( NX ) I / AL OG I 2 . ) ♦ . 5 
DO  2 IZ=1,NZ 
DO  21  I L= 1 * NX 
ESR( IL ) = EFR( IL  , IZ) 

21  ESMILI-EF  KIL.IZI 

CALL  FFT8A2(LNX,ESR,ESI ) 

DO  22  I L= I • NX 
EFR( IL, IZ)=ESR( IL) 

22  EFI ( IL, I Z I = ES I i IL  ) 

2 CONTINUE 
DSUM=0. 

L = 0 

DO  5 IZ=1,NZ 
SUM=0. 

DO  A I X= 1 , NX 

EFRlIX,IZ»  = SR2MErR( IX,IZ)+EFI( I X , I Z • » 

A SUM=SUMfEFR{ IX , IZ) **2 
SUM=SUM/FLOAT(NX) 

L*L  ♦ I 
SP(U=SUM 
IF(L.LT.A)  GOTO  5 
L * 0 

5 DSUMaDSUM+SUM 
OSUM=UEZ  *USUM 
WR I TE ( 6 , 91 ) DSUM 

90  FORMAT ( IH  , I 3 , AX , AE l A . 5 I 

91  FORMAT ( IHO , * Z INTEGRAL  OF  DC/C  SQUARED  = * , E I A . 5 , / ) 
IRSKP=INT(UX/DR*.l ) 

IZSKP=FL0AT(N)/FLUAT(NZ-1  >♦■.5 

RRSKP=I ./FLOAT { IkSKP) 

KZSKP=I ./FLOAT ( I ZSKP ) 

AZ ( I )=AZINZ)=YZ1 I ) =0. 

NZI=NZ-l 
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vSa>< 


100  CUNT INUE 

IF(SIG.EU.O)  RETURN 
IX=NUD(  1 IR-  1 »/  IRSKP  , NX  >4-1 
IXPl=IX*l  i I F l IX.EQ.NXIIXPl-1 
IF { IKSKP .UT . L ) GUTU  1 08 
00  104  I Z= 1 , NZ 
104  SST ( I Z ) =EFK ( I X , l Z i 

GU  TU  115 

108  JR=MJD ( IR-1 f IRSKP) 

U=FLUAT ( JK) *RRSKP 
V=l.-U 

OU  110  I Z = 1 » NZ 

110  SST  1 1 Z >=U*EFR( I XP 1 , 1 Z I ♦V*EFP|  I X, IZ  » 

115  CONTINUE 

DO  120  IZ*2,NZ1 

VZl lZ»=l./( 4.-YZI  I Z-  1 ) I 

TEMP=(  SST  C IZ*1)-2.*SSTJIZ)*SST(  IZ- 1 1 1 ♦RDBZS 
120  AZ  ( IZ  » = ( TEMP-AZI  IZ-1H  •YZII  L) 

DO  130  I B=2 ,NZ 1 
I Z = NZ ♦ 1 - I B 

130  AZ(  IZ)=AZl  IZI-YZl  IZI*AZ(  1 Z ♦ 1 I 
DO  140  I = l * NZ 1 
DO  140  J= 1 » I ZSKP 
IZ=IZSKP*( 1-1 )*J 
U=FLUAT(J)*RZSKP 
V=l.-U 

ssf  ( i z i =u*s  sti  ihi*v*sstiii- 

l DZ**2*U*V* ( AZ(  I I ♦ I 1 . ♦V)*AZ( I ♦ 1 » * « 1 . ♦ U » ) 

140  continue 

DU  150  1=1, N 
150  SSF ( I I=-FCTR*SSF(  I » 

RETURN 

END 


SUBROUTINE  TAPElWlIFLI 
COMMON/ AF AC TS/ AF ( 1 I ) 

COMMON/ NF  ACTS/BNFI 51 
COMMUN/MUNK/UM<  4> 

COMMON/ I WAVES/BIWI 1 1 I 
COMMON  / STEP/STE(4I 

COMMUN/ SAV/SAVRI I 00 , 4 ) , S A V II 100 , 4 » 
COMMON/ AC 0 US T/ UR ( 2 04  6 ) ,UI I 2048 > 
COMMON/VAR/R,lR,T,IT 
EQUIVALENCE  INSAV, STE ( 3» ) 

DIMENSION  BI2850I 
GO  TO  (100,200,3001  IFL 
100  CONTINUE 
L = 0 

00  1 1=1,11 
L = L*1 

1 BIL  >=AF( I ) 

00  2 1=1,5 
L = L ♦ 1 

2 b I L ) =BNF ( I i 
DO  3 1=1,4 
L = L*1 

3 B( L I = UM( I ) 

00  4 1=1,11 
L =L  ♦ 1 

4 BIL ! = BIW1  I ) 

DO  5 1=1,4 
L*L  ♦ 1 

5 B(L»=STE(  I I 

MKITEC 101  L.IBC I » , I = 1 ,L  » 

RETURN 

200  CONTINUE 
B( 1 l = R 
B(2»=T 
L = 2 

DU  31  1=1,1024 
L = L*l 

31  B( L I =UR ( I I 

DO  32  1=1,1024 
L = L*1 

32  Ql L » =U I ( I ) 

00  33  1=1,4 

OU  33  J=1 ,NSAV 
L«L*1 

B(  U=SAVR(  J,  I i 
L = L>  l 

33  B ( L I = SAV I ( J , I I 

WRITE ( 101  L, (Bill, 1=1, L) 

RETURN 

300  CONTINUE 
ENDFILE  10 
ENDfUE  10 
RETURN 
END 
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1 


SUbROUT I NE  UWS 

COMMON  /ACUUST/  UR ( 204 8 ), U II  2048  I 

COMMON  /PHASES/  SKU024» , CM  1 024  > , SMU024  ) ,CM(  1024) 

COMMON  /NFACTS/  N*LN»NR»NrLOSStNIT 

COMMON  /AFACrS/  FR E Q , SNORM , ZMAX , GW , FGW ( DR , FKR ( DZ , WL , P I , GD 

NM 1 =N- 1 

N2=N*2 

DU  10  J= 1 1 NM1 
UR ( J ♦ N I =-UK I N-J) 

10  UI ( J*N)=-UI IN-J1 
UR1N2»=0. 

UI (N2)=0. 

CALL  FFT8421LN*1 ,UR,U1  ) 

DO  20  J = 1 * NM 1 
JP1=J*1 

TEMP=CM  J )*URI  JP1 1-SKI  J )*U  I ( JPL) 

UI (JP1 )=-CK( J) *UI I JP1 I -SK ( J )*UR{ JP1 ) 

UR( JP1)=TEMP 
JN2=N2*1-J 

TEMP=CK( J )*UR( JN2 l-SKl J)*UI ( JN2 ) 

UI  ( JN2  ) =-CM  J ) *U  I ( JN2  ) -SK  ( J I *UR  ( JN2  1 
20  UR ( JN2 ) =T  EMP 
JR(1I=UI ( 1 )=0. 

U R ( N ♦ 1 I =U I ( N* 1 1=0. 

CALL  FFT8^2(LN*l ,UR,UI  ) 

00  30  J=  1 , N 

TEMP*CM( J)*UR< JI-SMI J ) *U I ( J) 

UI I J)=-CM( J )*UII J)-SM| J)*UR( J) 

30  URIJ)=TEMP 
RETURN 


5UBRGUT l Nt  ZMOOEJ 

COMMON  / I WAVES/  ENO , F [ NER , BZ , E L MAX ,NZ , AKM AX , NK , NMODE , 
1 SIG.TIW, JSTAR 

LARGE  EFR(256,257I ,EF  11256,2571 

COMMON  / EGNFN/  F ( 2 w 7 I , G< 2 5 7 I , QN { 2 57 I , QNS  ( 2 5 7 I 

COMMON  /VAR/  K,IR,T,IT 

REAL  TKl  257  I , ARC  256, 241,  Alt  256,241 

REAL  GS  V I 25b»  2 4 1 

REAL  KfKMAX»NORM»KCV 

DATA  PI, ERR  / 3 . L 4 1 5926 536 , I . E- l 0/ 

OAT  A ISTRT  / 0/ 

I F { ISTRT. NE.OI  GOTO  20 
I STRT=  I 
T P i =2  .*P l 

KMAX=AKM AX  % ZMAX=EZMAX 
DZ=ZMAX/(NZ-l  I $ DK=AKMAX/FLOATINKI 
NK2=NK*2 
UUM  -0 

DO  10  1=1, NZ 

z=dz*i  i-i  > 

TEMP=EN0**2*EXf><  -2 . *Z / BZ I -F I NER** 2 
TEMP=AMAXHO., TEMPI 
QN(NZ*l-l I * TEMP 
10  QNSINZU-I  l=SQRTI  TEMPI 

HN=ENO*BZ*l  l.-EXPl-ZMAX/aZI  I 
DO  12  I K= 1 , NK 
12  TK I I M =DK*( IK-  1 1 
20  CONTINUE 

IF!  IT.NE.il  GO  TO  26 
DU  25  MODE=l tNMUDE 
DO  24  I K*  l , NK2 
K=DK*  AMI.NO  I I K—  1 , NK  2 ♦ l - I K I 

P0H=-2.*PSIK,MU0£ I +ALOGI KANF I DUM l+l.E-1001 
AMP=SyRT I POW I 
ANG=TP l*RANF I DUM1 
AH  I IK,MODLI=AMP»COSIANGI 
24  AH  IK, MODE  I =AMP*S  1M  ANGI 
AR(NK*1, MODE  1=0. 

25  A I ( NK+ 1 , MODE  1 = 0. 

26  CONTINUE 

DO  27  IZ=l,NZ 
DO  27  I K= l , NK2 
EFR I IK, IZ 1=0. 

27  cF  II  IK, IZ 1 = 0. 

00  30  MODE  = 1 , NMODE 

GAMMA* I FLOAT ( MODE  I- .2  51 *P I /HN 

UG*OOG=0 . 

DO  30  I K= l , NK 
K=TKI 1KI 

GAMMA*  GAMMAS  OK  *1  UG*.5*DDG) 

IFIIT.NE.il  GAMMA=GSVI IK,MUDE1 
0G1=0G 

CALL  E lGENVL(GAMMA,K,NORM,ONORM,MODE,MERR, I ERR, 

•NZ ,OZ,HN,ERR I 

IFIMERR.EO.O.AND. IERR.EQ.OI  GU  TO  33 
KCT=500.*K/PI 

PRINT  302, MODE, KCY ,MERR, I ERR 


302  FORMAT  ( //40X,29Flfc  IGENVL  FAILED  TO  CONVERGE  AT,/ 

* 5 5X  » 5HMUDE  » I 2 » / <»5X  , 1 LHWAVENUMBEP  ,F  R.2/ 

♦ 45X , 7HMERR  = , I 2/ AbX , 7H I E RR  = ,I2> 

STUP 

33  KAT=NURM/UNORM 
DG=K*RAT/GAMMA 
DOG— Du- DG 1 

IF( IK.EQ.l)  DDG=UK*RAT /GAMMA 
IF(lT.EQ.l)  GSVI I K , MODE ) =GAMMA 
F I W = SORT  (FINER**2«-(  K/ GAMMA  ) **2 ) 

cw=cus(Fiw*n 

Sw=S IN ( F 1 w*  r ) 

IL=NK2*2-lK 
I F ( IK.EQ.l)  I L =NK ♦ 1 
AR I =AR ( IK, MODE) 

A I 1 = A I ( I K , MUDfc ) 

AR2  = AR ( I L ,MUDE ) 

A I 2 = AI ( I L , MODE ) 

DO  35  I 2 = 1 , NZ 
JZ=NZ*1-IZ 

EFR(IK,m=EFK(IK,IZ)  *F ( JZ ) * I AR 1 *CW- A I 1 *S W ) 

EF I ( IK , IZ )=EFI ( I K, I Z )*F< JZ ) »( AR1*SW*AI 1*CW) 

EFR ( IL  , I Z ) = EFR(  IL , I Z ) »F { JZ )*{ AR2*CW-AI 2*Srf > 

35  EF I ( IL, I Z >=EF I ( It , I Z ) ♦FI JZ ) *! AR2*SW*AI 2*CW) 

30  CONTINUE 
RETURN 
END 
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LOAD  MAP 


FI  REQUIRED  TO  LUAD  l<*7500 
HL  REQUIRED  TU  KUN  137000 

initial  transfer  to  user 


SLOCK  ASSIGNMENTS. 

SLOCK 

/NFACTS/ 

/AFACTS/ 

/VAR/ 

/ACOUST/ 

/STOP/ 

/PHASES/ 

/ BFUR/ 

/MUNK/ 

/ I rfAVES/ 

/SAV/ 

/STEP/ 

USER 
PKINTP 
T APE  I W 
DAYFOUL 
RNOMF 
/ EGNFN/ 

ZMQDEJ 

EIGENVL 

E 1 GENFN 

PS 

U*S 

I N I T 

PHASE  1 

PHASE2 

F F T 842 

ACGOEK 

ALNLUG 

ENDFIL 

EXP 

I BA l EX 

□UTPTB 

QUTPTC 

RANF 

RBAREX 

KUNSYS 

S 1 NCOS 

SORT 

STATUS 

TIMES 

GPUCIQ 

CPUSYS 

KODER 

Run i op 
SYSTM 
// 


A DDR  F S S 

LFNGTH 

100 

5 

105 

13 

120 

4 

124 

10000 

10124 

1 

10125 

10000 

20125 

6000 

2612  5 

4 

2o  1 3 1 

1 3 

26  144 

1440 

27604 

4 

27610 

1541 

31  351 

137 

31510 

5660 

37370 

170 

37560 

3135 

42715 

2004 

4472  1 

45154 

112075 

230 

1 12329^^^ 

236 

112563 

-44,1 

112714 

134" 

113050 

75 

113145 

54 

113221 

6161 

121402 

11242 

132644 

10 

132654 

1 10 

132764 

1 1 

132775 

53 

133050 

31 

133101 

31 

133132 

50 

133202 

24 

133226 

57 

133305 

303 

133610 

72 

133702 

42 

133744 

26 

133772 

37 

134031 

506 

134537 

53 

134612 

1235 

136047 

323 

136372 

322 

136714 

36 

FILE 


LGO 

LGU 

LGO 

LGU 

LGO 

LGO 
LGO 
LGO 
LGO 
4- GO 

LGU"- 

LGO 

LGO 

LGO 

SYSLIB 

SYSLIB 

SYSLIB 

SYSLIB 

SYSLIB 

SYSLIB 

SYSLIB 

SYSLIB 

SYSLIB 

SYSLIB 

SYSLIB 

SYSLIB 

SYSLIB 

SYSLIB 

SYSLIB 

SYSLIB 

SYSLIB 

SYSLIB 

SYSLIB 
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Appendix 

CUBIC  SPLINE  INTERPOLATION 


The  subroutine  RNDMF  uses  cubic  splines  to  interpolate  the  velocity 
proliles  from  the  coarse  mesh  on  which  the  random  internal  wave  field  is 
computed,  to  the  finer  mesh  on  which  the  acoustic  field  is  stored.  The 
theory  behind  this  is  as  follows. 

J 

/ 

We  want  an  interpolation  formula  for  t lie  function  v(z),  given  the 
mesh  points  (i  = l,2...,n+l)  and  function  values  v,  (1  = 1 ,2, . . . ,n+l) , 
bet 

AZi  = ZM-1  " Zi  ’ 1 = 1»2f..n  (A-la) 

AVi  = (Vi  + l 'Vi)/Azi  * i = 1,2 n . (A-lb) 

In  the  interval  z,  2 z i z , define  the  weight  functions 

i l + .l 

w = (z  - z^l/Az.  ( A-2 a ) 


w = 1 - w 

In  each  such  interval,  write 
the  mcsli  points  as  follows: 

v(z)  = wv  i-  wv 

i i + 1 

= wv  + wv 

i ill 


(A-2b) 


v(z)  as  a cubic  polynomial  continuous  at 


2 3 3 

+ ( A z. ) la, (w  - w)  + a (w  - w)] 
i i i+1 


- (Az.)  ww  Ta . (w+1 ) + a (w+l)| 

l i i + 1 


(A-3) 
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r 


rile  In  si  two  tonus  give  a linear  into  rpol  a t.ion  and  the  last  terms  give 
n second-order  correction  (cubic  in  z)  using  the  undetermined  coefficients 
n.  winch  arc  essentially  the  second  derivatives  at  the  mesh  points.  From 
Kq.  (A-3)  , we  find 


v'  (7.) 


A v + A 7 
1 i 


-a.  Uv'  - l)  + a ( 3 -l) 
1 i + 1 ' 


v"  (z)  = G 


a w + a w 
i i + 1 


(A-d ) 


(A-5 ) 


Equations  (A-3)  and  (A-5)  show  that  the  selected  way  of  writing  v(z)  makes 
both  v and  v continuous  functions  of  zr  (and  incidentally  Rq.  (A-5)  shows 
that  v"(z.)  =Gn.].  In  order  to  make  v' (z)  continuous,  we  match  values 

°f  V (ZiH)  u s i ng  Rq.  (A-d)  in  the  two  neighboring  intervals.  This  yields 
the  conditions 


Av , + A z . / a 1-  2a  \ = Av  - Az.  /2a  . \ 

1 Mi  Ml)  % ill  ^i  + i(2ji  + l H ai+2) 

which  can  be  put  in  the  form: 

Az  a -4-  2 (Az  h Az.  ) a . + Az  11 

\ i 1+1/  11I  i+1 

l‘»r  i = 1,2 n-1. 


1 . „ = A v -Av 

i+2  i+1  i 


(A -6) 


To  make  this  system  of  linear  equations  for  the  a determinate  we 
. i ’ 

choose  ai  = a^  - 0.  This  amounts  to  the  assumption  that  v(z)  is  linear 

in  z at  the  two  end  points  of  the  interpolation  region  (top  and  bottom  of 
the  ocean). 

Equation  (A-G)  is  a tri-diagonal  system  and  is  most  efficiently 
solved  numerically  using  the  "double  sweep"  method.  Look  for  a solution 
of  the  form 


a . = r . - s a 

1 i+l  i+1  ill 


i = 1 , . . . , n 


(A- 7) 
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Substituting  Eq.  (A-7)  into  (A-8)  gives  recursion  relations  for  r and 

i 

s as  foil ows : 
i 


a = — 
1 P . 


(fiv  - Av  \ - Az  r 
\ i i-1 ) i-1  i 


Az 


Pi  iJ  1 


(A-8) 


where  p.  = 2(Av  Av  ) - Az  s , i = 2 n 

1 l-l  i i-1  i 


Combining  Eqs.  (A-7)  and  (A-8)  we  obtain 


Az 


i+1  P, 


, i = 2. . . . ,n  , 


(A-9) 


i + 


Av  -Av  - Az  r 
i i-1  i-1  i 


/P  , i = 2 n 


(A-10) 


The  starting  values  of  these  right-sweeping  recursion  relations  are 
= s^  = 0 (as  follows  from  a ^ = 0).  The  full  procedure  is  therefore 
the  following.  Starting  with  r = s = 0,  compute  p from  Eq.  (A-8), 

A A A 

then  compute  s and  r from  Eqs.  (A-9)  and  (A-10),  and  again  use  Eq.  (A-8) 

A J 

to  compute  p , then  Eqs.  (A-9)  and  (A-10)  to  compute  s and  r , and  so  on 
3 4 4 

until  one  obtains  s and  r ; then  use  the  left-sweeping  recursion 

n+1  n+1  K 

relation  in  Eq.  (A-7)  starting  with  a = 0 to  compute  successively  a , 

n+ 1 n 

an-l a3'  a2  ‘ 

The  values  of  a^  computed  in  this  way  are  then  used  in  Eq , (A-3)  to 
compute  interpolated  values  of  v(z)  where  needed. 


A special  case  of  this  method  is  when  Az  = const  (say,  Az  = Az) 

i i 

for  all  i . Equations  (A-8),  (A-9),  and  (A-10)  can  then  be  simplified 
to: 

1 


i+1  4-s 


(A-ll ) 
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^ 01  * “ ' • . • i *'  • hf| nation  (A-7)  is  then  used  to  compute  n 

i n , . . . , 2 ; and  I>ci . (-A-.lt ) j^ives  v(z)  for  any  z. 


( A- 1 2 ) 
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FIGURE  1 DETERMINISTIC  SOUND-SPEED  PROFILE  AS  A FUNCTION  OF 
OCEAN  DEPTH.  The  value  of  c at  the  minimum  (zA  = 1300  m) 
•s  c = 1500  m/s. 
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FIGURE  2 SOUND-SPEED  PROFILES  DUE  TO  INTERNAL-WAVE 
MODES.  Realistic  internal-wave  spectra  have  significant 
intensities  for  horizontal  wave  number,  k,  less  than  about 
0.5  cycles/km.  Note  that  the  major  internal-wave  contri- 
butions to  sound-speed  fluctuations  occur  at  depths  less 
than  1 km. 
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FIGURE  3 


INTERNAL-WAVE  SPECTRUM  AS  A FUNCTION  OF  MODE 
NUMBER  j AND  HORIZONTAL  WAVE-NUMBER  k.  Although 
large  mode  numbers  contribute  verry  little  to  the  overall  spectrum, 
they  are  crucial  to  understanding  acoustic  effects,  since  their  vertical 
structure  allows  them  to  act  as  a scatterer  of  acoustic  energy  more 
readily  than  the  relatively  structureless  low  modes. 
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FIGURE  4 


TYPICAL  SOUND-SPEED  FLUCTUATION  PROFILE  INDUCED  BY  INTERNAL 
WAVES 
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FIGURE  5 ANALOGY  BETWEEN  (a)  AN  ACOUSTIC  RAY 
PROPAGATING  IN  A SOUND  CHANNEL 
CHARACTERIZED  BY  c (z),  AND  (b)  A 
SCHRODINGER-EQUATION  WAVE  PACKET 
IN  A POTENTIAL  WELL  V(z) 
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FIGURE  6 ANALOGY  BETWEEN  (a)  AN  ACOUSTIC-WAVE 
FUNCTION  TRAVELING  DOWN  THE  SOUND 
CHANNEL  WITHOUT  SPREADING  IN  DEPTH, 

AND  (b)  A SCHRODINGER-EQUATION  GROUND- 
STATE  WAVE  FUNCTION  IN  A POTENTIAL  WELL 
V(z).  Time  in  the  Schrodinger  situation  corresponds 
to  range  in  the  acoustic  case. 
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FIGURE  8 CODE  ORGANIZATION 
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